Discrete Feynman-Kac formulas for branching random walks 
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Branching random walks are key to the description of several physical and biological systems, such 
as neutron multiplication, genetics and population dynamics. For a broad class of such processes, 
in this Letter we derive the discrete Feynman-Kac equations for the probability and the moments 
of the number of visits ny of the walker to a given region V in the phase space. Feynman-Kac 
formulas for the residence times of Markovian processes are recovered in the diffusion limit. 



Consider the walk of a particle starting from a point- 
source and performing random displacements. At given 
times, the particle disappears by giving rise to k = 
0, 1, 2, ■ ■ • descendants with probability pk, the case pq 
corresponding to an absorption. Each descendant be- 
haves exactly as the mother particle, the overall path 
resulting in a branched structure, as shown in Fig. [T] 
Branching random walks lie at the heart of physical and 
biological modeling [U-IH , and are key to the description 
of neutron multiplication and nucleon cascades Q , disor- 
dered systems Q , evolution of biological populations @ , 
diffusion of reproducing bacteria 0, @], and mutation- 
propagation of genes [9[ , just to name a few. Detailed ac- 
counts of the huge research efforts devoted to this subject, 
ranging from the pioneering work by Galton and Watson 
on the extinction probability to the more recent develop- 
ments, can be found, e.g., in the monographs 
particular, even the simplest of these models, namely, a 
Brownian particle that at exponentially distributed times 
gives rise to two descendants, turns out to be highly non- 
trivial, and is still widely investigated, in view of its both 
practical and conceptual interest, especially in relation 
with such issues as front propagation and extreme statis- 
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A central question for random walks is to assess the so- 
journ time ty of the stochastic paths in a given portion V 
of the phase space Til- 19|: for branching Brownian mo- 



tion, this subject was first explored in connection with 
problems in mathematical genetics [i,0,Ei- The evo- 
lution of a number of physical and biological systems, 
such as neutrons or populations, is most often described 
in terms of discrete generations, so that one is quite nat- 
urally led to consider the number ny of visits to the do- 
main V, as illustrated in Fig. [U When V is large with 
respect to the typical displacement size, we can safely as- 
sume ny oc ty, which roughly amounts to assuming that 
the underlying walk can be approximated by a Brown- 
ian motion. However, it is well known that this simple 
proportionality breaks down when ny is small, and ap- 
plication of the diffusion approximation might therefore 



lead to inaccurate results [22j-|24|. 



In this Letter, we derive a discrete Feynman-Kac ap- 
proach to characterize the distribution P n (ny \tq) and the 
associated moments for a branching random walk that is 
observed up to the n-th generation, for arbitrary displace- 



ments, offspring distributions, and geometries. This pro- 
vides a general framework for dealing with a broad class 
of discrete-time random walks, which are ubiquitous in 
physics and biology. 

Feynman-Kac formulas. A single walker is initially 
isotropically emitted from ro, and undergoes a sequence 
of displacements and collisions. We assume that the 
displacements from r' to r, between any two collisions, 
are equally distributed, and obey the probability den- 
sity T(r';r). At each collision, the incident particle dis- 
appears, and k particles are emitted isotropically with 
probability pk ■ For the sake of simplicity, we furthermore 
assume that pk is spatially homogeneous. We formally 
define the number ny of visits to V as ny(n) = J^. V(ri), 
where V(ji) is the marker function of the region V, which 
takes the value 1 when the point r; G V, and vanishes 
elsewhere, and the sum is extended to all the points vis- 
ited by the source particle and its descendants up to the 
n-th generation. We assume here that the source point 
ro is not counted. Clearly, ny is a stochastic variable, 
depending on the realizations of the underlying process, 
and on the initial condition ro. The behavior of its dis- 
tribution, P n (ny |ro), is most easily described in terms of 
the associated generating function 
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FIG. 1: (Color online). Branching random walks. Left: the 
evolution of a path as a function of the generation number. 
Right: a path starting from ro and performing a number nv 
of visits to the region V . 
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This approach was first proposed by Kac, based on Feyn- 
man path integrals, for continuous-time Markov pro- 
cesses (see, e.g., [25|, f26|), and recently extended to 
non-Markovian walks [27| . The quantity i^ n (u|ro) is 
a stochastic functional (over the branching paths) that 
can be thought of as the discrete Laplace transform of 
P n (ny|ro). In order to derive an evolution equation for 
-F n (u|r ), we closely follow here the argument in [24|. We 
initially consider trajectories starting with a single parti- 
cle entering its first collision at ri. At r\, k particles are 
created, with probability pk ■ Random flights are Marko- 
vian at collision points, which allows splitting each of 
the k subsequent trajectories into a first jump, from ri 
to ri + Afc (the displacement obeying the jump length 
density T), and then a branching path from ri + A& to 
the positions held at the n-th generation. If k = 0, the 
trajectory ends at ri and there will be no further events 
contributing to ny. Hence, we have 

F n+1 («|n)=M- nri) + Pl u- v ^(F n (u\r 1 +A 1 )) + 

+ P2 u- v( ^(F n (u\r 1 + A' 2 )F n (u\ ri + A£)) + • • • ,(2) 

where expectation is taken with respect to the random 
displacements A, and u~ v ^ T1 ' can be singled out because 
it is not stochastic. The tilde is used to recall that we 
are considering trajectories starting with a single particle 
entering the first collision at ri . The terms in Eq. <j2j) can 
be understood as follows: the probability that k identical 
and indistinguishable particles (born at ri) give rise to 
ny collisions in V is given by the convolution product 
that the first makes n% collisions, the second n%, and 
the fc-th ny — n\ — n% — ■ ■ ■ . In the transformed space, 
this convolution becomes a simple product of generating 
functions. If we assume that the offspring particles are 
independent, the expectation of the products in Eq. (|5J) 
becomes the product of the expectations. We make then 
use of the discrete Dynkin's formula 



</(r 1+ A)) = J r*(r';ri)/(r')dr', 



(3) 



where / is any sufficiently well behaved function of a 
stochastic process, and T*(r';r) is the adjoint kernel as- 
sociated to T(r';r) Q. Intuitively, T* displaces the 
walker backward in time. We therefore obtain the dis- 
crete Feynman-Kac equation 



F n+ i(u|ri 



r*(r';ri)F„(u|r')dr' 



(4) 



where G[z] — X^fc^o Vkz k is the generating function of pk- 
Equation (jlj plays a central role, in that it relates the 
generating function F n of the number of visits ny to the 
generating function G of the offspring number k. Finally, 
by observing that the first collision coordinates ri obey 
the probability density T(r ;ri), it follows 



(5) 



F n (u\r a )= / F„(u|ri)T(r ;ri)dri 



together with the initial condition Fi(u\ti) = u~ v ( ri >. 
Knowledge of F n (u\r ) allows determining P n (ny | r ) : in- 
deed, F n (u\ro) is a polynomial in the variable u, the co- 
efficient of each power u~ l being P n (ny = i\ro). 

Moments formulas. A complementary tool for char- 
acterizing the distribution P n (ny\vo) is provided by the 
analysis of its moments. By construction, F n (u\vi) is the 
(rising) factorial moment generating function for trajec- 
tories entering their first collision at ri , which implies 



(n^Wri) = (-1)' 



d" 



du 7 



(6) 



= x(x + l)...(x + k— 1) being the rising factorial [28| . 
The tilde is again used to recall that the moments re- 
fer to trajectories starting from the first collision at ri. 
Combining Eqs. (Q| and © and using Faa di Bruno's for- 
mula [28| for the m-th derivative of the composite func- 
tion G[F n (r\)] yields the recursion property for m > 1 

(4 m3 ) n+1 ( ri ) = m V(r 1 )(h v n - 1) ) n+1 (r 1 ) + 



l V In 



(ri),- 



(m-j + l)\ 



)n(ri) , (7) 



where B m j[z-i,--- , ^ m -j+i] are the Bell's polynomi- 
als (28|, vj = (k(k— j+ 1)) are the falling facto- 
rial moments of the offspring number, and we have used 
rT*fr';ri)(n^)»(r')dr' = (n^) n (n). Bell poly nomi- 
als [31( commonly appear in connection with the com- 
binatorics of branched structures [28( : this might give a 
hint about their role in Eq. 0, which relates the mo- 
ments (hy^^n of the visit number to the moments Vj of 
the descendant number. Observe that (fiy ) n depends 
only on i>\, (h v ) n on V\ and V2, and so on. The recur- 
rence is initiated with the conditions (ny } n (r±) = 1, and 



Am) 



i(ri) = m\V(r±). Finally, the factorial moments 
)„(ro) for particles emitted at ro are obtained from 



(m) t 



n(ro) 



„(ri)T(r ;ri)dri 



(8) 



Stationary behavior. When trajectories are observed 



up to n 



we define 
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it turns out that the asymptotic moments 
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related to the stationary distribution of the particles [23] , 
provided that the limit exists. To see this, we introduce 
the collision density ^(r|ro) such that <I>(r|ro)<ir is the 
particle number at equilibrium in a volume of size dr 
around r, for a single particle emitted at r 0, [24| . It 
can be shown [29[ that the collision density satisfies the 
Boltzmann-like stationary integral transport equation 



f(r|r ) = ^ J T(r';r)*(r'|ro)dr'+r(r ; 



(9) 



where v — v\ is the average number of descendants per 
generation. Branching processes in the absence of spatial 
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FIG. 2: (Color online). Exponential flights on the interval 
[— R, R]: first factorial moment of ny as a function of n. Here 
Xo = 0, a = 1 and R — 1. Blue circles: ^ = 1.4 (J? c — 1.59); 
red stars: v = 1.2 (i? c ~ 2.57); green triangles: v = 1; cyan 
crosses: ^ = 0.8. Dashed lines: the stationary moments. 



FIG. 3: (Color online). Exponential flights on the interval 
[— R, R]: second factorial moment of nv as a function of n. 
Here xo = Q, a = 1 and R — 1. Blue circles: v — 1.4; red 
stars: v = 1.2; green triangles: v = 1; cyan crosses: ^ = 0.8. 
Dashed lines: the stationary moments. 



constraints are said to be subcritical for v < 1, supercrit- 
ical for v > 1, and critical for ^ = 1, respectively [H, 0, H[ . 
Equation can be understood as follows: the equilib- 
rium particle density at r for a source emitting at ro is 
given by the sum of all contributions having a collision at 
r', being multiplied and then transported to r, plus the 
contribution of the particles emitted from the source and 
never collided up to entering r. Now, from the proper- 
ties of the Bell's polynomials, the m-th order term can be 
singled out of the sum in Eq. ([7J , so that we can rewrite 
the stationary moment equation for (fiy 1 ^) as 

(n^Xrt) - v [ T*(r';ri)(^ m V)rfr' = (10 ) 



where <? m (ri) = mV(ri)(hy 1 )(ri) + b m (ri) acts as a 
source term, and 

m 

represents the contributions from Pk>2, and vanishes 
when po + Pi = E The integral equation (flU)) for 
{fly ){y\) can be explicitly solved (see [24|), and com- 
bined with ([5]) yields 

(4 m) )(r )= / *(r'|r ). 9m (r')dr', (12) 



starting with gi(r') = V(r'), which provides the desired 
relation between the stationary moments and the colli- 
sion density. When po + Pi = 1 one recovers the station- 
ary moment formula in 23] . 



An example. We illustrate this formalism on exponen- 
tial flights, a stochastic process that is widely adopted to 
model the propagation of neutrons and chemical or bio- 
logical species, among others (see, e.g., [2^, [23, Ho[ ) • In 
one dimension, the displacement kernel reads T(x')x) = 
a -i e -\x-x \/a ^ w ] lere <j j s some typical length scale. We 
assume that the source is io £ V. To fix the ideas, 
we take V as the interval [— R, R], the particles being 
lost upon crossing the outer boundaries. A natural ques- 
tion is whether the total number of visits will eventu- 
ally explose or level off to an asymptotic value under the 
combined effects of the branching mechanism and the 
spatial leakages. To our best knowledge, the full distri- 
bution P n (nv\^o) for this example is unfortunately not 
known. Actually, Eq. (UJ) can seldom be inverted, and 
the nonlinear integral equation resulting from taking the 
limit F = linin^oo F n hardly admits explicit solutions. 
Nonetheless, the moment Eqs. ([7]) and (fT2")) are amenable 
to exact formulas (at least for simple geometries and dis- 
placement kernels, such as in this case) and therefore 
greatly help in the analysis of branching random walks. 
Indeed, the essential features of a stochastic system are 
often captured by the first few moments. The resulting 
expressions are fairly cumbersome and will therefore not 
be reported here (hints on how to perform the calcula- 
tions can be found, e.g., in |29(). Instead, we display 
the first (Fig. [2]) and second factorial moment (Fig. [3]) 
as a function of the generation number n for some de- 
scendant number probability pk- Analytical results have 
been verified by Monte Carlo simulation. When v < 1 the 
moments of ny converge to an asymptotic value; when 
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v > 1, the moments may converge or diverge depending 
on whether the leakages from the boundaries compensate 
the growth of the particle number. This intuitively de- 
pends on the existence of a stationary solution ^(x\xq): 
we can then define a critical radius as the smallest R c 
such that for R > R c the average number of visits to V 
diverges when n — > +00. This quantity can be computed 
explicitly once (n^)(ro) is known, and for v > 1 reads 
R c = a esc -1 (y[v)l \/v — 1, independent of Xq. 

Diffusion limit. We conclude by examining the scal- 
ing limit of the discrete Feynman-Kac equations, which 
is obtained when ny is large, and at the same time the 
typical jump length e as well as the net average displace- 
ment \x are vanishing small. We set ty = nydt and 
t = ndt, where dt is some small time scale related to 
e by the usual diffusion scaling e 2 = 2Ddt and [i = vdt, 
D playing the role of a diffusion coefficient, and v of a 
velocity. By properly taking the limit of large ny and 
vanishing dt, ty converges to the residence time in V. 
It is natural to set pk = Xkdt, the quantity A& being 
a rate per unit of dt. Observe that when e and [i are 
small for any displacement kernel we have the Taylor ex- 
pansion / T(r'; r)/(r')dr' ~ /(r) - M 5 r /(r) + \e 2 d 2 J(v), 
and similarly / T*(r'; v Q )f{v')dv' ~ /(r ) + /ic\ /(r ) + I 11 
^e 2 d 2 a f(r ). It is expedient to introduce the quantity 
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Qt(u\r a ) = F t (e"|r ), which is the moment generating 
function of ty, i.e., 



<i?) t (r ) = (-1)' 



?t(«|r )| u =o, (13) 



when trajectories are observed up to time t. Under the 
previous hypotheses, combining Eqs. (j4]) and ([5]) and 
passing to the limit dt — > yields 



ggt 
dt 



= L* To Q t -uV{v )Q t + \G [Q t ] 



(14) 



where C* a = Dd 2 Q + vd ro — A, and A = J2k ^k- Finally, 
from Eq. (|13[) stems the recursion property 



^ (^) t (r )+mV(r )(C" 1 ) t (ro) + 



dt 



-\^VjB m ,j (t v ) t (r ), ■ ■ ■ ,{t 
3=1 



V I 



t(ro) .(15) 



As a particular case, when G[z] = 0.5 + 0.5z 2 and v = 
we recover the results in (2l| for critical binary branching 
Brownian motion without drift. 

Perspectives. The Feynman-Kac fomalism proposed 
in this Letter could be further generalized to take into 
account spatial dependences (T(r'; r) and/or pk may vary 
as a function of r) and anisotropics. A nice application 
of these results would be to infer the distribution pk on 
the basis of the number of countings ny recorded at a 
detector. 
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